Data from the test plate from IMR using “Universal Primers” from Parada et al, 2016. I followed the pipeline from Jesse McNichols and the Fuhrman lab protocols page.
Questions:
read=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/raw_seqs.csv',header=T,row.names=1)
data=data.frame(read[,7:8])
#barplot(data$Raw)
Same input for cDNA for all community samples (25ng into 20 uL reaction).
sub=subset(read, Frac_Comm=="Comm")
par(mar=c(5,5,0,1))
barplot(sub$Raw,las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,log='x',
xlab="Log10 No. Raw Reads")
Fractions had different RNA inputs, so different cDNA concentrations
sub=subset(read, Frac_Comm=="Frac")
par(mar=c(5,5,0,1))
barplot(sub$Raw,las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,log='x',
xlab="Log No. Raw Reads")
No real correlation between cDNA/RNA concentration and the number of reads per sample.
plot(log10(sub$Raw), sub$Frac.RNA, xlab='Log10 No. of Raw Reads', ylab='cDNA concentration', pch =21, bg='black')
Percent abundance of prokaryotic, eukaryotic, and cyanobacterial reads after sorting raw sequences based on hits to Silva and PR2 databases
On average, 1.8% of the total sequences were Eukaryotes
On average 0.46% of the raw reads were eukaryotes, so about 163/47,000 reads.
reads=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/pct_euk_bact.csv',header=T,row.names=1)
dim(reads)
## [1] 95 17
data1=data.frame(reads[,13:15])
sub=subset(reads, F_C=="Comm")
data1=data.frame(sub[,13:15])
par(mar=c(7,5,1,1),xpd=T)
barplot(as.matrix(t(data1)),las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,
xlab="% Abundance", col=c('navy','forestgreen','firebrick'),space=0)
box(which='plot')
legend(0,-10, legend=c("Prok", "Cyano", "Euks"), pch =22, pt.bg=c("firebrick", 'forestgreen','navy'),
bty='n', ncol=3)
On average, 3.4% of the fractionated sequences were eukaryotes, so about 806/24,000 reads.
Percentage of total raw reads that match to prokaryotic, eukarytoic, or cyanobacterial databases. The samples that don’t reach to 100 did not have matches to the databases, therefore we can assume they are probably chimera.s
sub=subset(reads, F_C=="Frac")
data1=data.frame(sub[,13:15])
barplot(as.matrix(t(data1)),las=2,names.arg = row.names(sub),horiz=T,cex.names=0.5,
xlab="% Abundance", col=c('navy','forestgreen','firebrick'),space=0)
box(which='plot')
legend(0,-10, legend=c("Prok", "Cyano", "Euks"), pch =22, pt.bg=c("firebrick", 'forestgreen','navy'),
bty='n', ncol=3)
Total of 110 unique ASVs.
library(phyloseq)
x<-read.csv(file='/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/feature-table.biom.csv',header=TRUE,row.names=1)
OTU = otu_table(x, taxa_are_rows=T)
taxa<-read.csv(file="/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/taxonomy.csv",header=TRUE,row.names=1)
t<-as.matrix(taxa)
tax2<-tax_table(t)
map<-import_qiime_sample_data("/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/map.txt")
#tree=read.tree("/Users/ahern/R/trophic_cascades/mc2/4Nov21/ASVs/tree.nwk")
phyo = phyloseq(OTU, tax2, map)
phyo
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 110 taxa and 84 samples ]
## sample_data() Sample Data: [ 84 samples by 10 sample variables ]
## tax_table() Taxonomy Table: [ 110 taxa by 9 taxonomic ranks ]
After QC - Average of 340 euk reads/sample + 3.4 ASVs per sample overall
library(ggplot2)
#phyo_abund=subset_samples(phyo, Comm_Frac=="Comm")
phyo_abund=transform_sample_counts(phyo, function(x) x/sum(x))
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class_18s.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = ASV)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "none")
d_rgn
phy=subset_taxa(phyo, Class!="Fungi")
#phyo_abund=subset_samples(phyo, Comm_Frac=="Comm")
phyo_abund=transform_sample_counts(phy, function(x) x/sum(x))
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class_18s.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
#sub=subset_samples(phyo_abund, Comm_Frac=="Frac")
#write.csv(otu_table(sub),'18S_nofungi.csv')
#write.csv(tax_table(sub),'18S_nofungi_tax.csv')
Note that the bar on the right is the heaviest of the fraction.
We can begin to see differences in abundance at the heavier vs. lighter fractions with even a few 18S samples
sub=subset_samples(phyo_abund, Comm_Frac=="Frac")
control <- subset_samples(sub,SampleID == "Ace.TB5.1" | SampleID == "Ace.TB5.4" | SampleID == "Ace.TB5.5" | SampleID == "Ace.TB5.15" | SampleID == "C12.TB5.7" | SampleID == "C12.TB5.11" | SampleID == "Met.TB5.3" | SampleID == "Met.TB5.9" | SampleID == "Met.TB5.10" | SampleID == "Met.TB5.12")
pd <- psmelt(control)
#colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class_18s.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, BD), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
Relatively the same, except my data does not have any unknowns, potentially due to under sequencing.
Sequencing Depth (seqs/sample):
Note: Opisthokonta is Fungi
#colors=randomcoloR::randomColor(110)
#write_rds(colors, 'colors_class.rds')
colors <- c("#bfd3e6","#fa9fb5", "#74c476", "#fc8d59", "#807dba", "#238443","#e31a1c", "#ec7014", "#969696")
read=read.csv(file='controlotu.csv',header=T,row.names=1)
colors <- c("#bfd3e6","#fa9fb5", "#74c476", "#fc8d59", "#807dba", "#238443","#e31a1c", "#ec7014", "#969696")
par(mar=c(15,5,1,10),xpd=T)
barplot(as.matrix(read), las=2, col = colors, space =0,
ylab="% Relative Abundance")
box(which='plot')
legend(8,.7, legend=row.names(read), pch=22,
pt.bg=colors)
Total of 1,632 ASVs, 1,608 ASVs if you get rid of Mitochondria and chloroplasts.
library(phyloseq)
x<-read.csv(file='/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/16S/feature-table.biom.csv',header=TRUE,row.names=1)
OTU = otu_table(x, taxa_are_rows=T)
taxa<-read.csv(file="/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/16S/taxonomy_16s.csv",header=TRUE,row.names=1)
t<-as.matrix(taxa)
tax2<-tax_table(t)
map<-import_qiime_sample_data("/Users/oliviaahern/Documents/MBL_WHOI/Trophic_Cascades/Experiment4/Sequences/IMR_Jan2023/18S/map.txt")
#tree=read.tree("/Users/ahern/R/trophic_cascades/mc2/4Nov21/ASVs/tree.nwk")
phyo = phyloseq(OTU, tax2, map)
phyo1 = subset_taxa(phyo, !Order==" Chloroplast")
phyo2 = subset_taxa(phyo1, !Family==" Mitochondria")
phyo=phyo2
agg=tax_glom(phyo2, taxrank="Class")
t=transform_sample_counts(agg, function(x) x/sum(x))
sub=subset_samples(t, Comm_Frac=="Comm")
phyo_abund=subset_samples(t, Comm_Frac=="Comm")
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(52)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class.rds')
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
Two samples did not amplify:
phyo_abund=subset_samples(t, Comm_Frac=="Frac")
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(52)
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, BD), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
Most Frequent Classes:
read1=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/fracs_class.csv',header=T,row.names=1)
data=data.frame(t(read1))
control=subset(data, Treatment=="Control")
meth=subset(data, Treatment=="Methanol")
acetate=subset(data,Treatment=="Acetate")
{
par(mar=c(5,5,1,1),mfrow=c(2,2))
plot(control$BD, control$Gammaproteobacteria,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),
xlab="Buoyant Density", ylab="% Abundance Gammaproteobacteria")
lines(meth$BD, meth$Gammaproteobacteria,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$Gammaproteobacteria,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
legend(1.76,.3, legend=c("Control", "Methanol", "Acetate"), text.col=c('gray70',"#a44f9a","#56ae6c"), bty='n')
plot(control$BD, control$Campylobacteria,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),xlab="Buoyant Density", ylab="% Abundance Campylobacteria")
lines(meth$BD, meth$Campylobacteria,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$Campylobacteria,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
plot(control$BD, control$Alphaproteobacteria,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),xlab="Buoyant Density", ylab="% Abundance Alphaproteobacteria")
lines(meth$BD, meth$Alphaproteobacteria,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$Alphaproteobacteria,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
plot(control$BD, control$Bacteroidia,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),xlab="Buoyant Density", ylab="% Abundance Bacteroidia")
lines(meth$BD, meth$Bacteroidia,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$Bacteroidia,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
}
abund=transform_sample_counts(phyo, function(x) x/sum(x))
phyo_abund=subset_samples(abund, Comm_Frac=="Comm")
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(1610)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class.rds')
colors_asv=c('gray70', colors)
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors_asv))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
abund=transform_sample_counts(phyo, function(x) x/sum(x))
phyo_abund=subset_samples(abund, Comm_Frac=="Frac")
pd <- psmelt(phyo_abund)
#colors=randomcoloR::randomColor(1610)
#write_rds(colors, 'colors_class.rds')
colors=readRDS('/Users/oliviaahern/Documents/GitHub/TestPlate/colors_class.rds')
colors_asv=c('gray70', colors)
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors_asv))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
Most Frequent ASVs
phyo_abund=transform_sample_counts(phyo, function(x) x/sum(x))
comm1=subset_samples(phyo_abund, Comm_Frac=="Comm")
s=subset_taxa(comm1, strain=="16e2e2408cedf8f9faae33004fb5eaa2" | strain=="4785f91e7f34f642ba5880ff1f2e7526" | strain=="ddb2bf017eeb122792d652cc207cb3b0" | strain=="639b8d71f6806556291f6114837ecc41" |
strain=="72043b4ee056432fd25f465446b2ff27" | strain =="d1c8f7b6346015d2d94b7db3e6250ffc")
pd <- psmelt(s)
colors_asv=c("#801900",
"#ffc77e",
"#525b00",
"#b8c636",
"#013c9e",
"#dc74e6")
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = strain)) +
scale_fill_manual(values=as.character(t(colors_asv))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn
Some facts to note
read1=read.csv(file='/Users/oliviaahern/Documents/GitHub/TestPlate/16S_asv_subset.csv',header=T,row.names=1)
control=subset(read1, Treatment=="Control")
meth=subset(read1, Treatment=="Methanol")
acetate=subset(read1,Treatment=="Acetate")
{
par(mar=c(5,5,1,1),mfrow=c(2,3))
plot(control$BD, control$X16e2e2408cedf8f9faae33004fb5eaa2,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),
xlab="Buoyant Density", ylab="% Abundance Aeromonas")
lines(meth$BD, meth$X16e2e2408cedf8f9faae33004fb5eaa2,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$X16e2e2408cedf8f9faae33004fb5eaa2,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
legend(1.76,.8, legend=c("Control", "Methanol", "Acetate"), text.col=c('gray70',"#a44f9a","#56ae6c"), bty='n')
plot(control$BD, control$X4785f91e7f34f642ba5880ff1f2e7526,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),
xlab="Buoyant Density", ylab="% Abundance Cellvibrio")
lines(meth$BD, meth$X4785f91e7f34f642ba5880ff1f2e7526,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$X4785f91e7f34f642ba5880ff1f2e7526,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
legend(1.76,.8, legend=c("Control", "Methanol", "Acetate"), text.col=c('gray70',"#a44f9a","#56ae6c"), bty='n')
plot(control$BD, control$ddb2bf017eeb122792d652cc207cb3b0,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),
xlab="Buoyant Density", ylab="% Abundance Chlorobi")
lines(meth$BD, meth$ddb2bf017eeb122792d652cc207cb3b0,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$ddb2bf017eeb122792d652cc207cb3b0,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
legend(1.76,.8, legend=c("Control", "Methanol", "Acetate"), text.col=c('gray70',"#a44f9a","#56ae6c"), bty='n')
plot(control$BD, control$X639b8d71f6806556291f6114837ecc41,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),
xlab="Buoyant Density", ylab="% Abundance Methylotenera 1")
lines(meth$BD, meth$X639b8d71f6806556291f6114837ecc41,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$X639b8d71f6806556291f6114837ecc41,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
legend(1.76,.8, legend=c("Control", "Methanol", "Acetate"), text.col=c('gray70',"#a44f9a","#56ae6c"), bty='n')
plot(control$BD, control$X72043b4ee056432fd25f465446b2ff27,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),
xlab="Buoyant Density", ylab="% Abundance Methylotenera 2")
lines(meth$BD, meth$X72043b4ee056432fd25f465446b2ff27,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$X72043b4ee056432fd25f465446b2ff27,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
legend(1.76,.8, legend=c("Control", "Methanol", "Acetate"), text.col=c('gray70',"#a44f9a","#56ae6c"), bty='n')
plot(control$BD, control$d1c8f7b6346015d2d94b7db3e6250ffc,
type='o', pch=21,col='gray70', bg='gray90', ylim=c(0,1), xlim=c(1.76, 1.83),
xlab="Buoyant Density", ylab="% Abundance Rhodobacter")
lines(meth$BD, meth$d1c8f7b6346015d2d94b7db3e6250ffc,
type="o", col='#a44f9a', pch=21,
bg="#a44f9a")
lines(acetate$BD, acetate$d1c8f7b6346015d2d94b7db3e6250ffc,
type="o", pch=21, bg="#56ae6c", col="#56ae6c")
legend(1.76,.8, legend=c("Control", "Methanol", "Acetate"), text.col=c('gray70',"#a44f9a","#56ae6c"), bty='n')
}
sub=subset_samples(phyo, Repeats=="Yes")
otu=otu_table(sub)
library(vegan)
spec=specnumber(t(otu))
par(mar=c(5,10,1,1))
barplot(spec, horiz=T, las=2, xlab='No. of Species')
There appears to be PCR bias in the nested PCR.
cycles=c(25,35,40,50,55,25,35,40,50,55,25,35,40,50,55)
par(mar=c(5,5,1,1))
plot(cycles, spec, pch=21, bg='black',
xlab="No. of PCR Cycles", ylab="No. of Species")
abund=transform_sample_counts(phyo, function(x) x/sum(x))
phyo_abund=subset_samples(abund, Comm_Frac=="Frac")
sub=subset_samples(phyo_abund, Repeats=="Yes")
pd <- psmelt(sub)
colors_asv=c('gray70', colors)
library(grid)
#palette(as.character(t(colors)))
d_rgn = ggplot(pd, aes(x = reorder(Sample, Timepoint), y = Abundance, fill = Class)) +
scale_fill_manual(values=as.character(t(colors_asv))) +
geom_bar(stat = "identity", width=0.97, color='black',
lwd=0.1) +
facet_grid(~Treatment, scale='free_x', space="free", shrink=TRUE) +
labs(x=" ", y = "Class Relative Abundance") +
theme_bw() +
theme(axis.text = element_text(color ='black',size=12, hjust=1,vjust=1),
axis.text.x=element_text(size=8, angle=90),
axis.title = element_text(color='black',face='bold',size=14),
legend.position = "bottom",
panel.grid=element_blank(),
strip.text.x = element_text(
size = 10, color = "black", face = "bold"),
legend.text = element_text(size=6),
legend.key.size = unit(0.25, 'cm')
#strip.background = element_rect(
# color="white", fill="white", size=1, linetype="solid"),
# panel.spacing = unit(0.05, "lines"))
)
d_rgn